rm(list = ls())
library(ggplot2)
library(tidyverse)
library(patchwork)
# 参数设置
kp <- 4
theta <- 0.36
gm <- 0.9
K0 <- c(0.8*kp^(1/(1-theta)),kp^(1/(1-theta)),1.2*kp^(1/(1-theta)))

# 变量初始化
lambda <- Y <- K <- matrix(NA,118,3) %>% as.data.frame()
lambda[1,] <- 1
K[1,] <- K0
Y[1,] <- lambda[1,]*K[1,]^theta*65^(1-theta)

# 动态追踪
for (i in 2:118) {
  lambda[i,] <- (1-gm)+gm*lambda[i-1,]+ runif(3,-0.02,0.02)
  K[i,] <- kp*lambda[i,]*K[i-1,]^theta
  Y[i,] <- lambda[i,]*K[i,]^theta*65^(1-theta)
}

# 画图
pK <- ggplot(K, aes(x = 1:118, y = V1)) + geom_line() + geom_line(aes(y = V2), linetype = 2) +
  geom_line(aes(y = V3), linetype = 3) + labs(x = '时间', y = 'K') + theme_bw()
pY <- ggplot(Y, aes(x = 1:118, y = V1)) + geom_line() + geom_line(aes(y = V2), linetype = 2) +
  geom_line(aes(y = V3), linetype = 3) + labs(x = '时间', y = 'Y') + theme_bw()
pK / pY
ggsave('../olg.png')
